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Abstract 

Micellization is the precipitation of lipids from aqueous solution into aggregates with a broad 
distribution of aggregation number. Three eras of micellization are characterized in a simple kinetic 
model of Becker-Doring type. The model asigns the same constant energy to the (k — 1) monomer- 
monomer bonds in a linear chain of k particles. The number of monomers decreases sharply and 
many clusters of small size are produced during the first era. During the second era, nucleii are 
increasing steadily in size until their distribution becomes a self-similar solution of the diffusion 
equation. Lastly, when the average size of the nucleii becomes comparable to its equilibrium value, 
a simple mean-field Fokker-Planck equation describes the final era until the equilibrium distribution 
is reached. 

PACS numbers: 82.70.Uv, 83.80.Qr, 05.40.-a, 05.20.Dd 
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I. INTRODUCTION 



Spontaneous self-assembly of small molecular aggregates in aqueous solutions forms as- 
sociation colloids or complex fluids ]!]. Depending on their mean aggregation number, 
molecular volume and critical hydrocarbon chain length, lipids can pack into spherical or 
cylindrical micelles. The surfaces of these structures are formed by the hydrophilic heads 
of the monomer molecules, whose hydrophobic tails lie inside the aggregate. Equilibrium 
thermodynamics shows that rod-like cylindrical aggregates have a polydisperse distribution 
of sizes (micellization), whereas the sizes of spherical aggregates grow indefinitely (phase seg- 
regation) 0. The latter process is similar to other examples of first order phase transitions 
such as condensation of liquid droplets from a supersaturated vapor, colloidal crystalliza- 
tion P| and the segregation by coarsening of binary alloys quenched into the miscibility gap 
[|. [5|, |6j . Understanding the kinetics of nucleation and growth beyond the determination of 
the steady-state nucleation rate is a task of great importance and not yet completely ac- 
complished. This is so despite the large literature on nucleation and growth 0, and several 
attempts at bridging the gap between nucleation and late-stage coarsening theories ||. 

In this paper, we study asymptotically a simple discrete model of micellization kinetics 
of Becker-Doring type |7], |8|, |9|, lCfl . Starting from an initial condition of pure monomers, 
we expect the system to evolve to the well known polydisperse equilibrium distribution [|TJ . 
However, the nonequilibrium evolution is interesting per se and because the methodology 
employed here may be applicable to the kinetics of phase segregation. We find that the 
approach to equilibrium occurs in three well defined stages or eras. Starting from the initial 
state of pure monomers, the number of monomers decreases sharply and many clusters of 
small size are produced during the first era. During the second era, aggregates are increas- 
ing steadily in size until their distribution becomes a self-similar solution of the diffusion 
equation. Lastly, when the average size of the nucleii becomes comparable to its equilibrium 
value, a simple mean-field Fokker-Planck equation describes the final era until the equilib- 
rium distribution is reached. Numerical solution of the model confirms all the theoretical 
predictions. 

The rest of the paper is as follows. In Section |IJ we review the equilibrium properties 
of self- assembling aggregates and introduce discrete kinetic models of Becker-Doring type 
to describe them. Depending on the binding energy of the aggregate with k monomers 
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(k cluster), micellization or phase segregation occurs. For rod-like aggregates, the bind- 
ing energy of a k cluster (relative to isolated monomers in solution) is (k — 1) times the 
monomer- monomer bond energy, and an equilibrium size distribution exists (micellization). 
For spherical aggregates, the binding energy includes a term proportional to the surface area 
of the aggregate and no equilibrium size distribution exists beyond a critical density. Then 
aggregates grow indefinitely and phase segregation occurs following the typical nucleation 
and growth kinetics. Section JT| presents a numerical simulation of micellization kinetics 
which clearly revels its three eras. The agenda of the asymptotic analysis is now clear, and 



is carried out in Section |V|. The last section contains our conclusions and suggestions for 
experiments. 

II. THERMODYNAMICS AND KINETIC MODELS 

The model presented here is nucleation in a lattice. There are systems, such as proteins 
aggregating in a cubic phase of lipid bilayers, for which a lattice formulation is physically 
correct. In this paper, the main reasons for a lattice model are clarity, and the expectation 
that the dilute limit of the lattice model (in which there are many more binding sites, 
M, than particles, N) should closely resemble crystallization from a dilute solution. The 
latter is a classical problem in the kinetic theory of first order phase transitions @. We 
shall now review the equilibrium statistical mechanics of aggregates, distinguishing between 
micellization and phase segregation, and then introduce the kinetic models we study. 

A. Equilibrium size distribution of aggregates 

Let us assume that we have p k > clusters with k particles (in short, k clusters) so that 

N 

N = Y,k Pk . (2.1) 

k=l 

Let efc be the energy of a A; cluster. The total energy of the lattice system is 

N N 

E = ^Pkek = N e 1 + ^p k (e k - hex), (2.2) 

k=l k=2 
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where we have used the particle conservation ( [2.1|) . Except for a constant Nei, the total 
energy is 

JV 

E = -Y,Pk£k, (2.3) 

k=2 

e k = kei-e k . (2.4) 

Now E is the total lattice energy measured with respect to a configuration in which all 
clusters are monomers, and e k is the binding energy of the k cluster (notice the sign conven- 
tion). We will obtain the equilibrium configuration by minimizing the free energy density 
with respect to the density of k clusters. To calculate the entropy, we proceed as follows. Let 
rij > be the occupation number of the site j, j = 1, . . . , M. The configuration space of the 
lattice consists of all M-tuples of occupation numbers, {ni, . . . , um}, with J2jLi n j = N and 
iV <C M. Clearly, there are many indistinguishable configurations that produce the same 
given set of numbers, pi, ■ ■ ■ ,Pn- Their number Q is given by the Bose-Einstein counting 
argument, 

_ Ml 

~ Pi\...p N l(M- Pl -...-p N )V (2 ' 5) 

and the entropy of the system is k B hxQ. In the appropriate thermodynamic limit, iV — > oo 
with fixed densities p = N/M (particles) and p k = Pk/M (k clusters), particle conservation 
becomes 

oo 
k=l 

and we can show that the entropy density is 

k]3 ( °° \ 

S = — lnfi ~ —k B ^2pkhipk + rlnrj , (2.7) 

oo 

r = l-X>, (2.8) 

k=l 

by using Stirling's formula. The free energy density, / = E/M — TS, can be written in 
terms of p and the densities of clusters having two or more particles by using its definition 



and Equations ( |2.3|) and (|2.6|) to (|2.8| ). The result is 



/ = - ^Pk£k + k B Tj2pk^Pk + k B Tr\nr, (2.9) 



k=2 k=l 
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where p\ = p — Y^T=2^Pk an d r = 1 — Y^=iPk- I n the dilute limit, 1 — r = J2kLiPk < 
SfeLi kpk = /)<!, and therefore r ~ 1, r lnr ~ — J2T=i Pk, an d Eq. (|2.9|) becomes 



/ = - E ^ + ^ T E ( ln /^ - 1), (2.10) 

fc=2 fc=l 

which corresponds to the Boltzmann counting. The equilibrium density of k clusters (k > 2) 
can be found by differentiating this equation with respect to pk and equating the result to 
zero. Taking into consideration that dpi/dp k = —k (k > 2), we obtain 

h = p " exp (SO (2,ll) 

(the positive sign in the argument of the exponential is due to our definition of the binding 
energies). Eq. ( |2.11| ) can be rewritten as 

9k = -e k + k B Tk\n(^j . (2.13) 

g k as a function of k can be interpreted as the activation energy of nucleation theory. The 
equilibrium density of monomers can be found by inserting Eq. Q2.ll ) into Eq. (|2.6|) and 
solving the resulting self-consistent equation for p\ in terms of the constant density p: 

Whether this self-consistent equation has a solution depends on the value of p and on the 
model we adopt for the binding energy of a A; cluster. Typical models are as follows. For 
rod-like aggregates, 

e k = (k - l)ak B T, (2.15) 
where ak B T is the monomer-monomer bonding energy [Q. For spherical aggregates, 

3 2 

£ k ~ (A; - l)ak B T - -ah? (2.16) 

for k ^> 1. Here a = 2^(Aixv 2 /3)^, where 7 and v = V/M are the interfacial free energy per 
unit area (surface tension) and the molecular volume, respectively. 

Inserting Eq. ( |2.15| ) in Eq. ( |2.14j ) and using J2T=i kz k = x ^Y!k=i xk = z /(l ~ x ) 2 > we 
obtain 



This equation has the unique solution 



_ 1 + 2pe« - 

Pl " 2^ ' (2 ' 18) 



with pi < e a for all values of the density p M. Notice that 



s = vi+w-i ( , 19) 

2^=1 Pk 2 

is the average cluster size in equilibrium. Notice that, for pe a ^> 1, (A;) ~ ^/pe" and 
p fc ~ e- a e~ k/{k) . 

For spherical aggregates, the self-consistency condition based on the approximation to e k 
in Eq. (|2~T^ ) is 

Pi £ k (p^- 1 exp = P- (2-20) 

Clearly this series converges provided pie a < 1 and it diverges if p\e a > 1. The critical 
monomer concentration p x = e _a is called critical micelle concentration (CMC) Below 
CMC eq. ( |2.20| ) can be solved for pi and the aggregates eventually form micelles with an 
equilibrium size distribution, whereas phase segregation and indefinite aggregate growth 
results if more monomers are added above the CMC. For > 1, the free energy ( ^.l^j) is 

2 

g k ~ ctksT + 3akz /2 — kip, with ip = kBT\n{pie a ). For ip > 0, g>& increases for small k, it 
has a maximum at the critical cluster size k c w (a/ip) 3 , and then it decays monotonically as 
/c further increases. 



B. Kinetic models 

Let us now formulate the kinetic theory of aggregation in these systems. As in the Becker- 
Doring kinetic theory, we shall assume that a k cluster can grow or decay by capturing or 
shedding one monomer at a time. Then 

Pk = Jk-i - jk = -D- jk, k>2, (2.21) 
jk = d k { e "b t p lPk - p k+1 \ , (2.22) 



or finally, 



3k = d k \[e ^ -1) P k- D +Pk \ , (2.23) 
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Here D + Sk = £k+i — £& — —D + gk + fesTln(l/pi) and j k is the net rate of creation of a 
k + 1 cluster from a /c cluster, given by the mass action law. We have made the detailed 
balance assumption to relate the kinetic coefficient for monomer aggregation to that of decay 
of a (k + 1) cluster, dk- Then pk given by Eq. ( |2.11| ) solves jk = 0. The kinetic model is 
described by a closed system of equations once we supplement Eqs. (pif), ( |2.21| ) and fl2.22| ) 



with expressions for the binding energy of a k cluster, Ek and for the kinetic coefficient of 
the decay reaction, dk. 

The simplest possible model for micellization within the Becker-Doring theory is obtained 
by setting Ek = {k — l)ak,BT and dk = 1 in Eq. (|2.22j ) for the creation rate of a {k + 1) 
cluster (rescaling of time can absorb a constant cluster decay rate dk = d; typical time 



scales describing aggregation kinetics range from microseconds to milliseconds [11, O, 1 



Equations ( |2.21 ) and ( 2.22Q then become the following discrete Smoluchowski equation: 



Pk + (e°Pi - 1) (pk ~ Pk-i) = Pk+i ~ 2pk + Pk-i, (2.24) 



to be solved together with the conservation condition ( |2.6| ), namely, Y^k^i^Pk — P- At 
t = 0, we assume that pk = p$ki- We shall consider the limit p ^> e~ a , in which the initial 

monomer concentration is much larger than the CMC. The parameters p and a are not 
really independent: if we rescale the cluster densities with p, so that 

Pk = prk, (2.25) 

and define a scaled time, 

r = e a pt= t - ) (2.26) 

the rescaled problem contains the single parameter e = (pe Q ) _1 1. Then Eqs. ( |2.24|) and 
( |2~6l) become 

dvu 

— ^ + (n - e) (r k - r fc _i) = e {r k+1 - 2r k + r fc _i), fe > 2 (2.27) 
dr 

oo 

1 = ^fcr fc , (2.28) 
fc=i 

to be solved with initial conditions 

n(0) = l, r 2 (0) = r 3 (0) = . . . = 0. (2.29) 
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Lastly, notice that we can straightforwardly derive two global identities from Eqs. (|2.27|) 
and ( ggg ): 

-, h rAri + r c ) + e (r x - r 2 - r c ) = 0, (2.30) 

ar 

dv 

— - + nr c + e (n - r c ) = 0. (2.31) 
ar 

Here r c is the total density of clusters 



oo 

r c = E r fc> ( 2 - 32 ) 



fc=i 

and, initially, r c (0) = 1. 

III. NUMERICAL RESULTS 

Numerical solution of the initial value problem given by Eqs. fl2.27j ) - ( |2.29| ) clearly ex- 
presses the phenomenology of micellization, and informs the singular perturbation analysis 
carried out in Section fV[ Figures |3J] to |3.5| illustrate the evolution of the size distribution 



for e = 4.54 x 10 4 (corresponding to a = 10 and p = 0.1). Figures |34](a), |3~2| , |3~3| and |3T4 



are histograms of as a function of k at different times, and Figure [3.5| records the time 
dependent behavior of the average cluster size, (k). 

Figure [3J](a) depicts an early stage of the kinetics. The sequences of small dots at each 
k record the values of at times between r = and r = 2, in increments of Ar = 0.2, and 
the larger dots joined by straight lines record the values of at r = 10. The direction of 
increasing time is generally clear: As indicated in Fig. |3.1| (b), the monomer concentration 
rapidly decreases to a small fraction of its initial value r\ — 1, so that the time orientation 
on the line k = 1 is downward. Many small clusters of sizes k, 2 < k < 5 are simultaneously 
created so the time orientation on the lines of these k is generally upward. Notice that p<i 
reaches a maximum and then decreases to a constant value, as can be seen in Fig. |3.1| (c). By 
the end of the initial stage at time r = 10, the creation of smaller clusters, with 2 < k < 5, 
has slowed down greatly relative to the initial spurt for times < r < 2. Furthermore, the 
number of clusters with more than 5 monomers is negligible. At r = 10, (k) m 2.69, much 
smaller than the equilibrium value (k) m y 'pe a = w 46.9. To determine the time scales 
appropriate for exploring the subsequent kinetics, it is highly instructive to plot the average 
cluster size (k) as a function of time, based on the numerical solution. Fig. ^T5] is a log-log 
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FIG. 3.1: (a) Scaled cluster size distribution pk/p as a function of k for < r < 10. At time r = 10, 
the values of pi/p, P2/P, etc. have been joined by straight lines as a guide for the eye. (b) Evolution 
of the scaled monomer concentration, p\j p. (c) Evolution of the scaled dimer concentration, pij P- 
Parameter values are a = 10 and p = 0.1. 



plot of {k)/e as a function of r. It reveals an initial rapid growth of (k) to a "plateau value" 
close to e, roughly located in the interval 10 < r < 100. In the subsequent growth after 
the plateau, large clusters with k 3> 1 eventually appear. Fig. [T5] indicates that by time 
t = 5 x 10 4 , k clusters having (k) w 10 are prevalent. 



Fig. 3.2 shows frames at times r = 20, 10 and 2 x 10 , thereby continuing those in 



Fig. The heavy dots correpond to r = 20, which is well inside the plateau phase. The 
histograms at r = 10 4 and 2 x 10 4 indicate the clear emergence of a continuum limit of the 
kinetics. 

In the time interval 2 x 10 4 < r < 5 x 10 5 , the log-log plot of (k)/e as a function of r 



in Fig. |3.5| is close to a straight line of slope 1/2. This strongly supports the existence of a 
self-similar stage of the kinetics. The line graphs in Fig. |3.4| depict {k) 2 rk as a function of 
x = kj (k) for the times r = 0.5 x 10 5 , 10 5 , 1.5 x 10 5 . They are nearly superimposed on top 
of each other. The heavy dots correspond to the plateau time r = 20, so the change in the 
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FIG. 3.2: Same as Fig. 3T(a), for the times r = 20, 10 4 and 2 x 10 4 . At the two later times, we 
have joined values of Pk/p corresponding to neighboring /c's by straight lines as a guide for the eye. 
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FIG. 3.3: Same as Fig. 3.1(a), starting at r = 2 x 10 5 . Snapshots of the size distribution have 
been taken at time intervals of r = 2 x 10 5 , until a time r = 16 x 10 5 . Then the last snapshot 
corresponds to r = 40 x 10 5 . 

distribution shape over the whole time span 20 < r < 1.5 x 10 5 is not very great. 

The self-similar stage is not the final chapter of the kinetics story either. By r = 10 6 , the 
linear dependence of ln((/c)/e) with lnr breaks down. In fact, at r = 10 6 , (k) « 31.1, which 
is comparable to the equilibrium value of 46.9 mentioned before. Evidently, there is a final 



stage of kinetics in which the size distribution asymptotes to its equilibrium form. Fig. 3.3 



is the final era of cluster aggregation, continued from Fig. |3.2| , in which snapshots of the size 
distribution are taken at r increments of 0.2 x 10 6 , from 0.2 x 10 6 to 4 x 10 6 . Convergence 
to an exponential distribution with (k) equal to the equilibrium value of 46.9 is clear. 
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FIG. 3.4: Approximate self-similar behavior of the size distribution at times r = 50, 000, 100,000 
and 150,000 (solid lines). Notice that (k} 2 rk is approximately the same function oik/(k) at different 
times. The dots correspond to r = 20. 
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FIG. 3.5: Evolution of the average cluster size, (k)/e, as a function of the scaled time r (thick 
solid line). The dotted line corresponds to the solution of the system ( 4.1 1| ) in Section [iy| below 



with an initial condition corresponding to the dot. The straight line of slope 1/2 corresponds to 
the self-similar continuum size distribution given by Eq. ( |4.19 ). 



IV. ASYMPTOTIC THEORY OF MICELLIZATION 

In this section, we shall interpret the numerical results shown in Section [TT1] by using 



singular perturbation methods; see Ref. [14] for a general description thereof. 
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A. Initial transient 



Initially, r\ (0) = 1 and there are no multiparticle aggregates. As we have seen in Section 
III], the numerical solution of the complete model shows that there is an initial transient 
stage during which dimers, trimers, etc. form at the expense of the monomers, and that 
rk ~ for sufficiently large k. Taking the e — > limit of Eqs. ( |2.30| ) and (|2.31|) yields the 
following planar dynamical system: 

dri , , 
-jj = -(n+rj, (4.1) 

Tr = ^ 
in the adaptive time s = J T p\dr. The general solution of the linear system ( |4.1| ) - ( |4.2| ) is 

T\ = (a — bs) e~ s , r c = be~ s , 

where a and b are arbitrary constants. Our initial condition yields a = b = 1, so that 

r 1 = (l-s)e" s , r c = e~% (4.4) 

and from Eq. (|4.3|) , 

rs e s 

ds. (4.5) 



lo 1 — s 

Clearly, r — > oo corresponds to s — > 1— . At s = 1, Eq. ( |4.4|) yields r± = 0, r c = e" 1 , which 
are the limiting values of the variables r\ and r c at the end of the initial stage. Eq. (|2.27 ) 
with e = becomes d(rke s )/ds = rk-ie s , which can be solved recursively to yield 



s k 1 s kS 



As t — ► oo, — > (A; — l)e _1 /A;!. Since r 6 (l) = 0.00255, after the initial transient stage there 
are insignificant numbers of aggregates with more than 5 monomers. In fact, the average 
aggregate cluster size is (k) = l/r c = e, whereas at equilibrium, (k) ~ y 'pe a ^> 1. We 
therefore conclude that there must be successive transients on time scales much larger than 
t = 0(e). 
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B. Intermediate transient 



Examination of the exact equation ( 2.27|) shows that when r± decreases to size 0(e) but 
r 2 , r 3 , ... are of order 1, all terms in its right hand side are 0(e). This suggests rescaling 
ri = eRi, so that p\ = e~ a R\, and using the original time t = er. Eq. ( |2.27| ) becomes 

^ = -(i?! - l)(r 2 - eR 1 ) + r 3 - 2r 2 + eR u (4.7) 
= -(R 1 - l)(r k - r fc _i) + r fe+1 - 2r fc + r k _ u k>2. (4.8) 



The global identities ( [2.30| ) and ( |2.31|) become 



( Rl -l)r c -r 2 + e + R\ + R 1 ] = 0, (4.9) 

dv 

-^ + (R 1 -l)r e + eR 1 = 0, (4.10) 

where now r c = eR\ + Y^h=2 r k ~ 11^=2 r k, as e ^ 0. In the limit e — > 0, Ri — 1 = r 2 /r c and 
Eq. ( [4.8| ) becomes 

dr k r 2 (r k -r k _ 1 ) 

— = hr fe+ i -2r k + r k -i, k>2. (4.11) 

at r c 

This is a closed system of equations for r 2 , r 3 , . . . , to be solved with the asymptotic values 
r k — (k — l)e~ 1 /A;! as initial conditions. It can be shown that the reduced versions of 
Eqs. (|4.10|) , r c = —(Ri — l)r c , and the conservation condition, J2kL2kr k = 1, are upheld 



automatically by the solution of Eqs. ( |4.11|) , so that they are redundant for this stage. 
The numerical solution of the reduced system of equations (|4.11 ) for r k , k > 2 closely 



approximates that of the full system of kinetic equations at this stage. It can be seen that 
more and more r k become different from zero as t increases and that r k — r k _i becomes 
small. This strongly suggests that r k can be approximated by a continuum limit for long 
times. To find the continuum limit, we set 

r k (t) ~ 5 a r(x, T), x = 5k,T = 5 b t. (4.12) 

Here 5-^0 fixes the scale of k = 0(1/5), so that x is fixed at some value of order 1. a and 
b are positive exponents to be determined. To find a, we use the conservation condition, 
ET =2 kr k = l: 

°° POO 

1 = S a ~ 2 Y(k5) r(k5, T) 6 ~ / x r{x, T) dx, 

k=2 Jo 

13 



provided a = 2. The limiting form of the particle conservation is thus 

xr(x,T) dx = 1. (4-13) 



o 



A similar calculation for the total number of clusters yields r c ~ 5 J °° r(x,T) dx, which 
suggests the definition 

rod 

r c ~5R c , R c = r(x,T)dx. (4.14) 
J o 



We now substitute Eq. ( }4.12|) in Eq. ( 4.11|) and use ( 4.14 ) instead of r c . The result is 



xb 9r S 2 r(26,T)[r(x,T)-r(x-6,T)}) 

5 — ~ — h r(x + 5,T)- 2r(x, T) + r(x - 5,T). 

ol o R c 

The right hand side of this expression is of order 0(S 2 ), so that the following distinguished 

limit is obtained if we set b = 2 and take 6 — > 0: 

dr(x,T) r(0,T) dr(x,T) d 2 r(x,T) 

dT R C (T) dx dx 2 ' 1 ' ' 

For k = 2, Eq. ( [4.1 1| ) and the scaling (|4 . 1 2|) with a = b = 2 imply that r(0, T) = 0. Therefore 

( |4.15| ) becomes the simple diffusion equation: 

dr d 2 r 

w = «?• < 416 » 

for x > 0, t > to be solved with the boundary condition r(0, T) = 0. 

The numerical solution of the discrete equations ( [4.1 1[ ) show that large aggregates do not 



emerge until t > 1. This suggests that the appropriate solution of Eq. ( [4.16| ) should be 



concentrated about x = as T — > 0+. That solution is proportional to the x derivative of 
the diffusion kernel, 



,2 , 

„2 



^-iH^r^^w (4l7) 

The numerical prefactor is chosen so that particle conservation, given by Eq. ( PI) , holds. 
It follows from Eq. ( |4.14|) that R c = (irT)~^. Hence the average aggregate size is 

<*> = *k = ■ (418) 

In terms of the original variables k, t and r&, the previous expressions are 

rt(t) ~ exp (4) ■ (4i9) 

(k) ~ v^rt, (4.20) 
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as t — > oo. These two equations yield 



(k) 2 r k 



irk 

W) 



cxp 



7T 

4 



k_ 

.<*>. 



(4.21] 



which resembles the behavior of the numerical solution of the full kinetic model as indicated 
in Fig. |3]4]. Notice that the average cluster size (k) corresponding to the solution of Eqs. 
Q4.1H ) (dotted line in Fig. p.5|) approaches the value Q4.20|) (straight line of slope 1/2 in Fig. 



C. Equilibrium transient 



The large time limit of Eq. (|4.19 ) does not match the equilibrium size distribution, which 
is r k ~ ee~ k ^ in the same scaled units; see Section [TT]. Thus the limit given by Eq. Q4.19 ) 



is expected to break down when it predicts an average (k) of the order of the equilibrium 
length l/y/e. According to Eq. (|4.20|) , this occurs at a time \/i = 0(e~2), i.e., t = 0(e _1 ). 
In this third and final transient towards equilibrium, we set: 

r k (t) = er(x,t), x = ^fek } T = et. (4.22) 



This is the same scaling as in Eq. (|4.12 ) with a = b = 2 and 5 = y/e, and therefore we use 



here the same notation for the variables. With this scaling, the scaled particle conservation 
is 

oo oo 

1 = ^^kr k = ea k r(x, T), 

k=l k=l 

and the limit e — > yields 

/ xr(x,T)dx=l. (4.23) 
Jo 

Similarly, 

i r°° i 
r c ~e2 / r(x,T)dx = R c . (4.24) 
Jo 

The scaled version of the global identity ( |2.32| ) is 

R c (R 1 -l) + e l iR l + e^ = 0. (4.25) 
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Here r\ = eR\ = er(e^,T). It follows from Eq. (4.25) that 



€1 



R 1 -l = -— + 0(e). 

K c 



(4.26) 



The scaled kinetic equation ( |2.27| ) is 



e 3 ^ = -e 2 (i?i - l)[r(x,T) - r(x - e*, T)] + e 2 [r(x + e*, T) — 2r(x, T) + r(x — e^, T)]. 

We now substitute Eq. ( |4.26| ) in this expression, divide it by e 3 and take the limit e — > 0. 
The result is 



dr 



1 dr d 2 r 
+ 



(4.27) 



dT R C (T) dx dx 2 ' 
In these units, the average aggregate length is (x) = l/R c , and Eq. Q4.271 ) can be rewritten 



as 



dr , dr d 2 r 
dT = ^ X 'dx + dx 2 ' 



(4.28) 



to be solved with the boundary condition 



r(0,T) = l, 



(4.29) 



which follows from Eq. ( [4.26| ) with e — ► 0. It can be straightforwardly checked that 



-iL J °° xr(x,T) dx = 0, and therefore / °° xr(x,T) dx = 1, provided r(x, 0) satisfies this 



particle conservation condition. 
We now have to show two things: 



1. As T — > 0+, the solution of Eqs. fl4.28|) and ( [4.29|) is asymptotic |14| to the right hand 
side of Eq. fl4.17|) , the self-similar limiting solution of the intermediate transient stage. 



2. The solution of Eqs. (|4.28| ) and fl4.29|) tends to the equilibrium size distribution as 
T — > oo. 

Then the size distribution of the equilibration transient as T — > 0+ matches the long time 
limit of the previous intermediate stage, and tends towards equilibrium as T — >• oo. This 
completes the description of the dynamics of the aggregate size distribution. 
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1. Matching with the intermediate transient stage 



We represent r(x,T) as 



r(x,T) = ±h((,T), ( = ^=. (4.30) 



With prefactor 1/T, the particle conservation equation ( |2.(j| ) and the total cluster density 
adopt the invariant forms 

/ (h((,T)d(=l, (4.31) 
Jo 

roc I roc h (J 1 ) 

Rc(T)=j Q r(x,T)dx = -^= j o h(C,T)dC=^J-. (4.32) 

Then 



<*> = m- (4 - 33) 

Inserting this equation together with Eq. ( |4.30|) in Eq. (|4.28| ), we obtain 

d 2 h , 18h _ fdh (h\ 

to be solved with the boundary condition indicated by Eqs. Q4.29|) and ( |4.30| ): 

h{0,T)=T. (4.35) 

Asymptotic similarity as T —>■ means that h((, T) in Eq. ( |4.3U| ) has a limit H(() as T — > 0. 
The limit equations obtained from Eqs. ( (4.31[ ), ( f4.34[ ) and ( |4.35| ) are 



d 2 H TT 1BH 

W + H+ 2 C 1( =0 m C> °' 

H(0) = 0, 



roo 

/ c^(C)dC=i- 

JO 



The unique solution of these equations is H(() = (e ^ 2 / 4 /(2-\/7r), which is the right hand 
side of Eq. ( p7|) . 



Trend towards equilibrium 



The stationary solution of Eq. ( [4.28 ) with the condition (|4.29|) is r e = e x( - x \ and the 



particle conservation condition gives (x) 2 = 1, so that (x) = 1. Then the stationary solution 
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of Eq. (|4.28|) is r e = e x , which is the sought equilibrium solution. To show that r(x,T) 
r e (x) as T — > oo, we define the associated free energy 

r 

To 



-r + r In 



<ix — 1, 



r 



(4.36) 
(4.37) 



and show that it is a Lyapunov functional for Eq. (|4.28|) . Notice that J °° r Inr^dx = 

— Jq°° xt dx = —1, and therefore f[r] is the usual free energy, f[r] = J °°(r lnr — r) dx. 
Firstly, the standard inequality xlnx > x — 1 for positive x = r/r yields / > 

— J °° e~ x dx — 1 = —2, and therefore / is bounded below. Notice that f[r ] = —2 at 
equilibrium. 

Secondly, time differentiation of Eq. ( [4.371 ) yields 



dT 



00 dr ( r 
o 9T ln U Ur 



If we now substitute Eq. (|4.28|) , integrate by parts, and use r(0, T) = r (0) = 1 and r dx 
l/(x), we obtain 

2 



df_ 
dT 



(x) 



00 1 / dr \ 2 
o r \dx I 



(x) 



oo r°° I ( dr 
rdx - — 
o jo r \ ox 



(4.38) 



The right hand side of this equation is less or equal than zero because of the Cauchy-Schwarz 
inequality: 



r(0,Tf 



00 dr \ 
o ax / 



<9r 



9x 



\ yoo r°° 1 f dr\ , 
dx ) < rdx - — ax. 
/ Jo Jo r \dx I 



Therefore, we have proven that the free energy is a Lyapunov functional. We can rewrite 
Eq. ( |4.38D in an equivalent form by defining f = exp[— x(x)], and using the identities: 

(x) = (x) 2 / rdx = [ r ( — ^— - \ dx, 
Jo Jo 



f°° dr 
(x) = — (x) / — dx 
Jo dx 



to obtain 



df _ 


POO 

— I r 


" a 




dT 


Jo 


dx 




This equation shows that r - 


-> r as 


T - 


-> oo. 



<9x 

91nr 9 lnf 

9x <9x 



rfx < 0. 



dx, 



(4.39) 



/ °° xr dx = 1 yields (x) 2 = 1, and therefore r = e 
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FIG. 4.1: Comparison of the approximation ( 4.40|) (dashed line) to the numerical solution of the 



full kinetic model (solid line) for four different times r: (a) 100,000, (b) 500,000, (c) 1,000,000, and 
(d) 3,000,000. Notice that the agreement improves as the equilibrium distribution is approached. 

3. Approximation of the size distribution function by matched asymptotic expansions 

An uniformly valid approximation to the size distribution function can be easily formed 
from: (i) rjj. (r), given by Eqs. ( |4.5| ) and (|4.6|), (ii) r^\t), which solves the approximate 
system of Equations ( f4.11|) and r c = X)j£L 2 r k with the initial conditions rfc(0) = (k— l)e~ l /k\, 
and (iii) r(x, T), which solves the nonlinear Fokker-Planck equation ( f4.28| ) with the condition 
( |4.29| ) and it matches ( |4.17| ) as T — > 0+. The result is 



r ( r f \r) = r£V) + 4V) + er(V&, eV) - - * exp (-f-) . (4.40) 

k\e 2A/7T(er)2 \ 4er/ 



Figures 4J. compare the distribution function given by Eq. ft4.40|) to the numerical solution 



of the complete model equations in times corresponding to the end of the intermediate stage 
and the beginning of the equilibration stage. At these times, = {k — l)/(k\e). We 
observe a good agreement between approximate and numerical solutions, which improves as 
the time elapses and the equilibrium distribution is approached. 
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V. CONCLUSIONS 



On the basis of a simple kinetic model and starting from the initial state of pure 
monomers, we have shown that the process of micellization of rod-like aggregates at high 
CMC occurs in three separated stages or eras. In the first era, many clusters of small size are 
produced while the number of monomers decreases sharply. During the second era, aggre- 
gates are increasing steadily in size and their distribution approaches a self-similar solution 
of the diffusion equation. Before the continuum limit can be realized, the average size of the 
nucleii becomes comparable to its equilibrium value, and a simple mean-field Fokker-Planck 
equation describes the final era until the equilibrium distribution is reached. A continuum 
size distribution does not describe micellization until the third era has started: during the 
first two eras the effects of discreteness dominate the dynamics. 

In order to validate our theory by an experiment, it would be important to measure 



the average cluster size as a function of time as in Fig. [3.5| : the multiscale behavior is 
more clearly seen in this figure. To determine the time scale, we need a measure of the 
cluster diffusion coefficient, d, which was set equal to 1 in Section [D]. A convenient relation 
could be Eq. ( [4.20| ), which in dimensional units is (k) « yd/nt. In case the self-similar size 
distribution is not reached during the intermediate phase, another way to determine d is 
to study the equilibration era and compare the experimentally obtained size distribution 
with the numerical solution of the model. At equilibrium, (k) 2 rs pe a , and this relation 
determines the dimensionless binding energy a. 
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